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Large-eddy simulation of flow around 
an airfoil on a structured mesh 

By Hans-Jakob Kaltenbach AND Haecheon Choi 

1. Motivation and objectives 

The diversity of flow characteristics encountered in a flow over an airfoil near 
maximum lift taxes the presently available statistical turbulence models. This work 
describes our first attempt to apply the technique of large-eddy simulation to a 
flow of aeronautical interest. The challenge for this simulation comes from the high 
Reynolds number of the flow as well as the variety of flow regimes encountered, 
including a thin laminar boundary layer at the nose, transition, boundary layer 
growth under adverse pressure gradient, incipient separation near the trailing edge, 
and merging of two shear layers at the trailing edge. 

The flow configuration chosen is a NACA 4412 airfoil near maximum lift. The 
corresponding angle of attack was determined independently by Wadcock (1987) 
and Hastings & Williams (1984, 1987) to be close to 12°. The simulation matches 
the chord Reynolds number U^c/v = 1.64 x 10 6 of Wadcock’s experiment. 

2. Accomplishments 

2.1 Numerical method and SGS model 

The numerical method for solving the unsteady, incompressible Navier-Stokes 
equations is described in Choi et al (1993). Second-order spatial central differences 
on a staggered mesh axe combined with a semi-implicit time integration scheme. 
Formulation of the problem in terms of contravariant velocity components, weighted 
with the Jacobian, in conjunction with the staggered variable configuration leads to 
discretized equations that can be solved with the classical splitting approach. The 
resulting pressure Poisson equation is solved using FFT for the spanwise (periodic) 
direction and iterative methods for the remaining two-dimensional problems. The 
computational cost is about equally distributed between computation of the right- 
hand side and solving the Poisson equation at every substep of a third order Runge 
Kutta time integration. 

The implementation of the dynamic subgrid-scale model (Germano et al 1991) 
with least-square contraction (Lilly 1992) uses the spanwise homogeneity of the flow 
to obtain a model coefficient that is a function of streamwise and wall- normal coor- 
dinate only. We found that the dynamic procedure occasionally renders unrealistic 
negative coefficients in regions where the flow is laminar such as at the nose or in 
the potential flow region. In these regions, the negative values are the result of the 
dynamic model becoming ill-conditioned and have no physical significance. In the 
present simulations we prevented any form of backscatter by constraining the model 
coefficient to be always positive. 
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Table 1. Spacing along upper surface of airfoil. The last two columns show cell 
aspect ratio and ratio of boundary layer thickness to domain width for case A. 


On the present mesh, the CFL limit of 1.5 results in an average timestep of 
2 x 10 -4 c/{7oo. About 80 CPU-seconds on a Cray-C90 are needed to advance the 
solution over one timestep on a mesh of 638 x 79 x 48 = 2.4 x 10 6 cells. Therefore, 
simulation of one time unit c/Uoo requires 90 CPU-hours. In order to obtain smooth 
statistics the results have to be averaged over several time units. 

2.2 Computational domain and mesh layout 

The computational domain is a C-mesh with the outer boundary about three 
chord lengths away from the surface. At the outer boundary we specify the freestream 
velocity Uoo • As a consequence, the vertical velocity component (in a coordinate 
system aligned with the chord at 0° angle of attack) will be zero at the outer bound- 
ary. Therefore, the chosen configuration resembles more the flow around an airfoil 
inside of a wind tunnel with parallel walls than an airfoil in free air. Jansen (1995) 
has shown that, even with the walls located much closer, the presence of wind tun- 
nel walls mainly affects the flow in the nose region by increasing the suction peak. 
The pressure distribution in the rear part and the size of the backflow zone, how- 
ever, are only weakly dependent on whether the wind tunnel walls are included or 
not. A no-slip condition is enforced at the airfoil surface, and we use a convective 
(radiative) boundary condition at the outflow plane. 

Results from two simulations will be presented. The two cases differ only with 
respect to the spanwise domain width which is 0.05c in case A and 0.025c in case 
B. The spanwise spacing A 2 is the same with 48 cells in case A and 24 cells in 
case B, respectively. Main criterion for the choice of the spanwise domain size 
is the ratio of boundary layer thickness to domain width, which is tabulated in 
Table 1. As a consequence of the rapid growth of the boundary layer thickness on 
the suction side, this ratio, which is initially sufficiently small to capture several 
structures in the spanwise direction, exceeds one near the trailing edge. It is likely 
that the development of flow structures in the outer part of the boundary layer will 
be affected by the limited domain size. Comparison of cases A and B gives insight 
in the sensitivity of the simulation with respect to this parameter. 

The design of an adequate mesh involves several aspects. The most energetic 
eddies of the boundary layer have to be resolved. More or less general criteria 
have been developed for the mesh spacing in the case of wall bounded shear flows 
under zero pressure gradient. However, these criteria depend on the numerical 
method employed (Lund et al, 1995). Cabot (1994) found that for LES of turbulent 
channel flow based on second-order finite differences a spacing of Ax+ = 60 and 
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FIGURE 1 . Time series of spanwise velocity fluctuation between stations x/c — 0.24 
(bottom) and x/c = 0.98 (top) at about 5% of the local boundary layer height. 
Individual curves are separated by a vertical offset of 0.3 with the corresponding 
zero-lines located at 1.2, 1.5, ... 3.9. 

A = 15-20 is needed to adequately resolve the near wall structures. 

Little is known about the minimum spacing requirements for boundary layers 
which are close to separation. The mesh size in terms of wall units probably becomes 
less relevant in this case. About half of the 640 streamwise points were distributed 
over the upper surface, which guaranteed that the streamwise spacing was between 
1/3 and 1/5 of the local boundary layer thickness for most of the upper surface, see 
Table 1. The streamwise spacing varies considerably along the surface due to the 
boundary layer growth. Near the trailing edge, the grid was refined in x in order 
to resolve the merging of the two shear layers. No attempt was made to resolve the 
turbulence on the lower side of the airfoil. Spacings in terms of wall units based 
on the local skin friction as given in Wad cock’s experiment are given in Table 1. It 
is evident that the spacing in the present simulation is considerably coarser than 
what has been found to be necessary for channel flow simulations. However, as 
the boundary layer develops along the surface, the resolution criteria become less 
restrictive so that the flow in the rear part is much better resolved than in the front 
section. 

In the wall-normal direction we used a hyperbolic mesh generator (Chan, 1993) 
to distribute 79 layers of cells. The first line away from the wall was at about 
y + = 1, and over most of the surface there were between 20 and 30 points inside 
the boundary layer. 

2.S Difficulties arising from the high Reynolds number 

Centered difference schemes suffer from the emergence of grid-to-grid oscillations 
(2-A-waves, wiggles) when used for high Reynolds number simulations. Usually, 
the viscosity provided by the subgrid-scale model is sufficient to dampen these grid- 
to-grid oscillations. Several sources for 2-A-waves have been identified in the past 
(Gresho, 1981). They include high cell Peclet numbers in conjunction with large 
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Figure 2. Mean velocity profiles, normalized by £/«>, along upper surface. Sym- 
bols: case A , B , measurements by Hastings o and Wadcock + . 

streamwise gradients of the advected variable. This situation is typically encoun- 
tered near the nose and the trailing edge of the airfoil. Other sources are the outflow 
boundary (an artificial boundary layer is generated in the streamwise direction) and 
mesh stretching. As shown by Cain & Bush (1994), waves propagating into an in- 
creasingly coarse (fine) mesh are amplified (dampened) in a centered scheme. In 
our simulation we find that strong 2- A- waves appear near the nose and near the 
trailing edge. The wiggles appear almost exclusively in the streamwise coordinate 
direction. Part of these waves travel with different phase speed and cancellation 
occurs. However, other parts are steady and accumulate in time. These stand- 
ing waves contaminate the potential flow region after long integration times. It is 
difficult to assess to what degree the solution is contaminated by the presence of 
2-A-waves. On a staggered mesh, velocity components are averaged in order to 
obtain fluxes at cell faces. This averaging on a scale of the mesh cell can sometimes 
completely hide the 2- A- wave. For example, the convective term d(uu)/dx in the 
streamwise momentum balance is evaluated as 


1 

( U.+1 + u, V 

( m + m-i\ 2 

Ax 

V 2 / 

V 2 ) . 


The finite difference expression renders the same value independently whether an 
oscillatory part in the i-direction u, = (-l)'u a with zero mean and arbitrary am- 
plitude u a is added or not. Similarly, if a 2-A-wave in the i-direction is present 
in the v velocity component, it will not appear in the discrete approximation for 
d(uv)/dy. However, it will contaminate the term d(uv)/dx. Time averaged fields 
of velocity components show 2-A-waves in the potential flow region, but the pres- 
sure field is virtually free of wiggles. This indicates that the presence of 2-A-waves 
in the potential flow region may be tolerated to a certain degree since wiggle free 
streamlines in accordance with the pressure field can be reconstructed. 
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FIGURE 3. Pressure distribution around the airfoil. Symbols: LES , Wad- 

cock o . 

The strongest effect of 2-A-waves comes from the associated limitations for the 
computational timestep. Large amplitude wiggles in the wall normal velocity com- 
ponent in conjunction with rather fine wallnormal spacing cause high CFL numbers 
near the nose. The resulting timestep limitations are so severe that the simulation 
can not be carried out at an affordable cost. We therefore resorted to an ad hoc 
modification of the numerical scheme. In a small region near the nose (less than 
2% of the chord) we applied a 1:2: 1-filter in the streamwise and spanwise direction 
which efficiently eliminates all 2-A-waves. Filtering is equivalent to adding a di- 
rection dependent diffusion term to the equations. Justification for this procedure 
comes from the fact that the flow near the nose is laminar and filtering on a scale 
of the grid cell does not affect the flow physics. Additionally, the boundary layer 
in the experiments was tripped at a location around x/c = 0.02, thereby fixing the 
region of laminar-turbulent transition. We find that the flow spontaneously transi- 
tions as soon as the filter ends. In this sense, we control the location of transition 
by setting the streamwise extent of the region where the solution is filtered. The 
filter extended about 40 layers away from the wall and faded to zero over another 
15 layers. Unfortunately, this procedure changes the potential flow significantly. 
Because the mesh cells are rather large in the outer part of the domain, filtering on 
the grid scale is no longer negligible on the scale where the potential flow changes 
near the nose. Future simulations can easily avoid this problem by limiting the filter 
to the vicinity of the surface, i.e. it should end near the boundary layer edge. No 
attempt was made to dampen 2-A-waves in the trailing edge region where the flow 
is fully turbulent. Any filtering there would probably affect the flow physics. 

2.4 Simulation results and discussion 

Figure 1 shows time series of the spanwise velocity fluctuation w recorded at 
several stations along the upper surface of the airfoil. We observe a shift in the 
frequency which corresponds to the most energetic motions towards lower values as 
the recording station moves closer to the trailing edge. This is consistent with the 
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FIGURE 4. Boundary layer thickness 6 , displacement thickness 6 *, momentum 
thickness 0 and shape factor H = 6* / 9 along the upper surface of the airfoil. Sym- 
bols: LES, • Hastings, x Wadcock. 


increase of an inertial timescale (ratio of the boundary layer thickness to the edge 
velocity) as the boundary layer grows under the influence of the adverse pressure 
gradient. It becomes evident that the solution has to be sampled over several time 
units c/Uoo in order to obtain representative turbulence statistics for the rear part 
of the airfoil. 

Statistics were obtained by averaging the instantaneous flow fields in the spanwise 
homogeneous direction and in time over more than 2c/Uoo* Profiles of the mean 
velocity in a surface normal coordinate system are shown in Fig. 2. At the first two 
stations, the edge velocity is about 12% smaller than measured by Hastings. As 
mentioned earlier, this is a side effect from the filter which was applied in the nose 
region in order to eliminate 2-A-waves. Since filtering was limited to a region close 
to the surface, simulated and measured mean flow agree much better for distances 
greater than y/c = 0.06. Although a better match between simulated and measured 
edge velocity is desirable (and can easily be obtained by further reducing the dis- 
tance from the surface over which the filter is applied), we don’t expect turbulence 
statistics to be significantly affected. One reason is the observation that the simu- 
lated adverse pressure gradient matches the measured one over most of the upper 
surface, see Fig. 3. Filtering affects mainly the magnitude of the suction peak and is 
partially responsible for the offset in the simulated pressure distribution. Addition- 
ally, since wind tunnel walls were not properly considered in this simulation, the 
pressure distribution near the nose will deviate from the measured one, see Jansen 
(1995). The goal of the present study is to predict the boundary layer growth and 





Structured mesh LES of NACA^IZ 


57 



x/c = 0.40 0.59 0.66 078 0.92 


FIGURE 5. RMS of velocity fluctuations u' ,v* and w ', normalized with Uqq. 

Symbols: case A , B , u' o and v* □ from Hastings. In terms of relative 

magnitude, the three curves for each simulation are v f , w * , and u ' respectively. 


the amount of separation near the trailing edge. Accurate prediction of the suction 
peak is of secondary interest. 

Displacement and momentum thickness from the simulation lie in between the 
measurements of Hastings &; Williams (1984) and Wadcock (1987) upstream of 
x/c = 0.4, see Fig. 4. The experimental values differ by up to 40% as a result of dif- 
ferences in boundary layer tripping and Reynolds number. However, the measured 
shape factor H « 1.55 is similar in both experiments in the region x/c = 0.2... 0.4. 
Contrary to the experiment, H drops gradually in the simulation in the region 
x/c = 0.2... 0.4 and reaches values as low as 1.4. 

Since both experiments measure similar boundary layer growth and flow retarda- 
tion near the trailing edge, the flow development does not seem to be very sensitive 
with respect to the exact values of 6* and 6 of the turbulent boundary that develops 
behind the transition strip. Although the thickness of the simulated boundary layer 
is close to the measured ones in the front part of the airfoil, the underprediction 
of the shape factor in the simulation and the initially opposite trend (decline as 
opposed to a growth) indicates insufficiencies in the simulated boundary layer for a 
considerable part of the upper surface. This is not surprising since the resolution is 
so coarse that the near wall structures can hardly be resolved properly. Examina- 
tion of instantaneous flow fields close to the surface reveals a very streaky structure 
with typical spacings in the order of a few mesh cells. Similarly, spanwise two point 
correlations show a zero-crossing within 2-3 spanwise grid points for all near-wall 
locations upstream of x/c = 0.5, see Fig. 6. This indicates that the simulation has 
marginal resolution near the wall. Further evidence comes from the comparison of 
the present case with an earlier simulation which was a factor of 2 coarser in the 
streamwise direction and a factor of 1.5 coarser in the spanwise direction. The flow 
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FIGURE 6. Spanwise two-point correlations R uu ( ), R vv ( o ) and R ww 

( ) versus distance zjc for three stations along the upper surface. Bottom 

figures correspond to the near wall region, top figures to y/6 & 0.3. The y coordinate 
of the lower figures is shifted, i.e. y = —1.5 corresponds to a zero correlation. 


retardation and the boundary layer growth was significantly improved on the finer 
mesh. Therefore, further grid refinement and, subsequently, a better prediction 
of the boundary layer in the front region might lead to better agreement between 
simulation and measurements over the entire airfoil. 

The shear stress provided by the SGS model is an indicator for the role of the 
SGS model. The maximum contribution is about 15% of the resolved stress uv and 
is found in the front part of the airfoil where the resolution is coarse. Near the 
trailing edge, the SGS stress is negligible compared to the resolved Reynolds shear 
stress. The ratio of SGS eddy viscosity to molecular viscosity is about 20, which 
emphasizes the important role of the model for the kinetic energy budget. 

RMS values of the velocity fluctuations are shown in Fig. 5. Agreement between 
simulation and experiment is reasonable in the middle section of the airfoil. In a 
characteristic manner for an adverse pressure gradient boundary layer, the location 
of maximum rms values (and Reynolds shear stress) moves towards the outer part of 
the boundary layer. Also, the anisotropy of the fluctuations in the outer part of the 
boundary layer is greatly reduced. Substantial differences between simulation and 
experiment are indicated by the large discrepancy in simulated and measured rms 
values (and shear stress) near the trailing edge. It is unclear whether this mismatch 
is a local effect or rather a result of differences in (spatial) flow history between 
experiment and simulation. 

Results from cases A and B, which only differ with respect to the spanwise domain 
size, are surprisingly similar. Two-point correlations from the outer part of the 
boundary layer of case A do not drop to zero within half the spanwise width for 
locations downstream of x/c = 0.6, see Fig. 6. This means that the large scales 
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of motion are affected by the presence of artificial periodic boundaries. Since the 
limitations are much more severe in case B as compared to A, one would expect 
that both cases deviate in the rear part. Presently, it is not clear why the simulation 
is rather insensitive with regard to the domain width. Kaltenbach (1994) made a 
similar observation for a flow in a diffuser where the aspect ratio of the outlet duct 
was smaller than 0.5. Doubling the aspect ratio had only a small effect on the flow 
evolution. The cost for case B is about half that of case A. Further studies on the 
effect of grid refinement would be much cheaper if the domain width of case B turns 
out to be sufficient. 

3. Conclusions and future goals 

Wall resolving LES of flow around an airfoil has been demonstrated to be feasi- 
ble with present computers and standard numerical schemes for LES. Qualitatively, 
the simulation captures typical features of separating flows such as boundary layer 
retardation and drastic increase in Reynolds stresses. This demonstrates the capa- 
bility of the LES concept to deal with flows in complex configurations of immediate 
technical interest. However, the resolution provided was probably too coarse to 
adequately simulate the boundary layer in the first half of the airfoil. Although the 
resolution might have been adequate for the rear part, the overall agreement with 
measurements with respect to prediction of backflow is not satisfactory. History 
effects might play a role, and further studies should attempt to match better the 
integral boundary layer parameters of the experiment at an early station. Because 
of conservation properties, the use of centered difference schemes is very desirable 
in the context of LES. However, the emergence of 2- A- waves is a serious problem 
for the present high Reynolds number flow and needs further consideration, for ex- 
ample, usage of explicit filters as explored by Lund & Kaltenbach in this volume. 
Comparison of two cases with different domain width did not show significant sen- 
sitivity with respect to this parameter in the range considered. Future simulations 
should consider the effect of wind tunnel (top and bottom) walls by a corresponding 
modification of domain size and boundary conditions. 
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